Load libraries
#library(GmAMisc)
library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.0 ──
## ✓ ggplot2 3.3.2 ✓ purrr 0.3.4
## ✓ tibble 3.0.4 ✓ dplyr 1.0.2
## ✓ tidyr 1.1.3 ✓ stringr 1.4.0
## ✓ readr 1.4.0 ✓ forcats 0.5.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
library(readxl)
library(ggplot2)
#library(dplyr)
library(ggthemes)
library(RColorBrewer)
library(plotly)
##
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
##
## last_plot
## The following object is masked from 'package:stats':
##
## filter
## The following object is masked from 'package:graphics':
##
## layout
library(car)
## Loading required package: carData
##
## Attaching package: 'car'
## The following object is masked from 'package:dplyr':
##
## recode
## The following object is masked from 'package:purrr':
##
## some
library(ordinal)
##
## Attaching package: 'ordinal'
## The following object is masked from 'package:plotly':
##
## slice
## The following object is masked from 'package:dplyr':
##
## slice
library(emmeans)
library(multcompView)
library(broom.mixed)
## Registered S3 method overwritten by 'broom.mixed':
## method from
## tidy.gamlss broom
library(forcats)
#library(plotly)
library(boot)
##
## Attaching package: 'boot'
## The following object is masked from 'package:car':
##
## logit
library(egg)
## Loading required package: gridExtra
##
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
##
## combine
Bootstrapping functions
# functions for bootstrapping 95% confidence intervals around the mean
# get mean and 95% CI of values x via bootstrapping
boot_ci <- function(x, perms=5000, bca=F) {
get_mean <- function(x, d) {
return(mean(x[d]))
}
x <- as.vector(na.omit(x))
mean <- mean(x)
if(bca){
boot <- boot.ci(boot(data=x,
statistic=get_mean,
R=perms,
parallel = "multicore",
ncpus = 4),
type="bca")
low <- boot$bca[1,4]
high <- boot$bca[1,5]
}else{
boot <- boot.ci(boot(data=x,
statistic=get_mean,
R=perms,
parallel = "multicore",
ncpus = 4),
type="perc")
low <- boot$perc[1,4]
high <- boot$perc[1,5]
}
c(low=low,mean=mean,high=high, N=round(length(x)))
}
# get mean and 95% CI via bootstrapping of values y within grouping variable x
boot_ci2 <- function(d=d, y=d$y, x=d$x, perms=5000, bca=F){
df <- data.frame(effect=unique(x))
df$low <- NA
df$mean <- NA
df$high <- NA
df$n.obs <- NA
for (i in 1:nrow(df)) {
ys <- y[which(x==df$effect[i])] #pulls out all the e.g.location prefs
if (length(ys)>1 & var(ys)>0 ){
b <- boot_ci(y[which(x==df$effect[i])], perms=perms, bca=bca) #resamples with replacement the mean and ci for e.g. AJs pref for location, gives back the low mean and high
df$low[i] <- b[1]
df$mean[i] <- b[2]
df$high[i] <- b[3]
df$n.obs[i] <- b[4]
}else{
df$low[i] <- min(ys)
df$mean[i] <- mean(ys)
df$high[i] <- max(ys)
df$n.obs[i] <- length(ys)
}
}
df
}
# y = response score, x= response within group (e.g. bat spp)
#use relevant working directory
# gleaner_universe_r_7_19_19 <- read_excel("~/Dropbox/GLEANER UNIVERSE STUDY/data sheets/gleaner_universe_r_7.19.19.xlsx",
# na = "NA")
g_u <- read.csv("https://raw.githubusercontent.com/maydixon/katydids/master/gleaner_universe_r_08.15.22_github.csv", sep=",", header=TRUE, stringsAsFactors = F )
str(g_u)
## 'data.frame': 1933 obs. of 28 variables:
## $ Species : chr "Trin_nice" "Lonc_auri" "Micr_hirs" "Phyl_hast" ...
## $ Species_full : chr "Trinycteris nicefori" "Lonchorina aurita" "Micronycteris hirsuta" "Phyllostomus hastatus" ...
## $ ID : chr "Naranjo" "Lonc_auri_01.08.18_1" "Ghandi" "King Lear" ...
## $ Sex : chr "M" "M" "M" "M" ...
## $ Date : chr "4/8/19" "1/8/19" "11/8/18" "2/1/19" ...
## $ Set : chr "2" "3" "2" "?" ...
## $ Attempt : int 1 1 1 1 1 2 1 1 1 1 ...
## $ Playback : chr "H. frenatus" "Q. gigas" "Beetle flight" "P. pustulosus" ...
## $ Number : chr "20" "16" "17r" "2" ...
## $ Time : chr "0:00" "22:48" "23:09" "19:11" ...
## $ Response : int 0 2 1 0 0 0 1 1 0 0 ...
## $ Notes : chr "0" "1" " " " " ...
## $ Reviewed : chr "" "" "" "" ...
## $ Exclude : chr "" "" "" "" ...
## $ call_duration : chr "2.13" "15.735" "0.490975" "0.3785" ...
## $ peak_freq : num 3788 1900 22525 858 858 ...
## $ min_freq : num 1890 700 525 812 812 ...
## $ max_freq : num 11148 3900 58775 2935 2935 ...
## $ bandw : num 9255 3100 58200 2120 2120 ...
## $ entropy : num 0.392 0.179 0.491 0.148 0.148 ...
## $ peaks : num 7.25 1 7.75 1 1 1 1 1 1 1 ...
## $ pulses_call : num 10.8 66.5 NA 2 2 ...
## $ calls_callgroup : num 1 1 NA 1 1 1.5 1 1.5 1 2 ...
## $ sum_pulse_duration_callgroup: num 0.256 15.735 NA 0.384 0.384 ...
## $ duration_callgroup : num 2.13 15.735 NA 0.379 0.379 ...
## $ call_rate_used : num 54 15.97 NA 1.95 1.95 ...
## $ duty_cycle : num 0.0394 0.9854 NA 0.1944 0.1944 ...
## $ file_time_with_sound : num 0.425 58.15 43.65 20.475 20.475 ...
Check and clean data
Code “call types” like “katydid” and “cicada”
#Recode a variable with the categories
g_u<- g_u %>% mutate(Call_type = if_else(Playback %in% c("A. spatulata","C. wheeleri","D. gigliotosi" ,"T. subfalcata","V. zederstedi"), "Katydid",
if_else(Playback %in% c( "Acla sp.","Aclodes sp.", "Anaxipha sp."), "Cricket",
if_else(Playback %in% c("D. diastema", "D. ebraccatus", "H. fleishmanni", "P. pustulosus" ,"S. sila"), "Frog",
if_else(Playback %in% c("Q. gigas", "Z. smaragdina"), "Cicada",
if_else(Playback %in% c( "H. frenatus","T. rapicauda" ), "Gecko",
if_else(Playback %in% c("Beetle flight", "Frog hopping" , "Mouse rustles", "Moth sounds" ),"Movement cue",
if_else(Playback %in% c("White noise"), "White noise","Other"
))))))))
#order factor
g_u$Call_type<- factor(g_u$Call_type, levels = c( "Cricket" , "Katydid", "Cicada" , "Frog" ,"Gecko", "Movement cue", "White noise" ), ordered = TRUE )
unique(g_u$Call_type)
## [1] Gecko Cicada Movement cue Frog Katydid
## [6] Cricket White noise
## 7 Levels: Cricket < Katydid < Cicada < Frog < Gecko < ... < White noise
Note Moth sounds were only played to a few bats, should be removed for some analyses
0 - no response
1 - ears twitch in time with call
2 - orients toward sound when playing
3- flies towards speaker
4 - hovers within .5 m
5- lands on speaker
Recode response with 0, 1, 4, 5, -> Make a column that codes
whether bats: 0 - (did nothing (0),
1- ear twitched or came closer(1-3),
4- approached within a half meter (4),
5- or landed (5).
(collapses 1-3, these distinctions may not be important/ espcalatory
across spp.)
Also make a binary column that just codes whether or not a bat
landed. ->
or 0, 1, 2-3,4,5. With the idea that orient and comes closer can be
conflated.
#Recode a variable with 0, 1, 4, 5
#g_u<- g_u %>% mutate(Response_ApproachLand = if_else(Response %in% c("0","1","2" ,"3"), "0", if_else(Playback %in% c("4"), "4", if_else(Response %in%c("5"),"5", "other"))))
g_u <- g_u %>% mutate(Response_ApproachLand = case_when(
Response == 0 ~ 0,
Response == 1 ~ 1,
Response == 2 ~ 1,
Response == 3 ~ 1,
Response == 4 ~ 4,
Response == 5 ~ 5))
g_u <- g_u %>% mutate(Response_Land = case_when(
Response == 0 ~ 0,
Response == 1 ~ 0,
Response == 2 ~ 0,
Response == 3 ~ 0,
Response == 4 ~ 0,
Response == 5 ~ 1))
Palette for coloring by Call_type
# Create a custom color scale (using RColorBrewer)
library(RColorBrewer)
CallCatColors <- brewer.pal(7,"Set1")
# Palette rearranged
#CallCatColors <- c("#984EA3", "#4DAF4A", "#377EB8", "#FF7F00", "#E41A1C", "#FFFF33", "#A6A6A6")
# Alternative palette (OG powerpoint scheme)
CallCatColors <- c("#2A769E", "#4D6123" , "#1E846A" , "#D0892C" , "#973721", "#573AA2" , "#A3A3A3" )
#Alternative palette 2 (modified powerpoint scheme)
#CallCatColors <- c("#FDE84C", "#4D6123", "#19657E", "#EE964B", "#8C2635", "#6D3A47", "#A6A6A6")
# View palette
pie(rep(1,7), col=CallCatColors )
names(CallCatColors) <- levels(g_u$Call_type)
colScale_CallCat <- scale_fill_manual(name = "Call type", values = CallCatColors)
colScale_color_CallCat <- scale_color_manual(name = "Call type", values = CallCatColors)
#and then add the color scale onto the plot as needed
#summarize # of individuals per species
n_species <-g_u %>% group_by(Species_full) %>%
summarise(N_total = n_distinct(ID, na.rm=TRUE), N_Female = length(unique(ID[Sex == 'F'])), N_Male = length(unique(ID[Sex == 'M']))) %>% arrange(desc(N_total))
## `summarise()` ungrouping output (override with `.groups` argument)
n_species
## # A tibble: 15 x 4
## Species_full N_total N_Female N_Male
## <chr> <int> <int> <int>
## 1 Trachops cirrhosus 13 5 8
## 2 Lophostoma silvicolum 12 2 10
## 3 Artibeus jamaicensis 10 1 9
## 4 Micronycteris microtis 10 4 5
## 5 Phyllostomus hastatus 10 0 10
## 6 Micronycteris hirsuta 7 2 5
## 7 Phyllostomus discolor 7 1 6
## 8 Gardnerycteris crenulatum 6 3 3
## 9 Lonchorina aurita 5 1 4
## 10 Lampronycteris brachiotis 3 1 1
## 11 Trinycteris nicefori 3 0 3
## 12 Lophostoma brasiliense 2 1 1
## 13 Tonatia saurophila 2 0 2
## 14 Glyphonycteris daviesi 1 0 1
## 15 Micronycteris schmidtorum 1 0 1
#sum total n bats
total_n <- sum(n_species$N_total)
total_n
## [1] 92
#write.csv(n_species, "n_species.csv")
bootstrapped means and CIs for playback for each bat spp (for faceted plots)
# # end goal: bootsums for all bats in all playbacks
# #but blanks for n<7 in high and low columns
# #put blanks in n<2 for mean column
#
# # n per species
# n <- g_u %>% group_by(Species_full) %>% summarize(n = length(unique(ID)))
#
# #bootstrap playback values within each bat spp
# Boot_sums <- g_u %>%
# filter(Call_type !="NA") %>%
# mutate(effect= paste(Species_full, Playback, sep="_")) %>%
#
# boot_ci2(d = ., y = .$Response, x = .$effect, perms = 5000, bca = F) %>%
# mutate(Species_Playback = effect) %>% #keep combined column in new name
# separate(effect, into=c('Species_full', 'Playback'), sep="_", remove=F) #re separate combined column
#
# # add n spp
# Boot_sums <- left_join(Boot_sums, n) # n per spp.
#
#
# # delete bootstrapped values below minimum sample size for plotting (could also do by total bat n, not observations per each playback)
# N_min_mean <- 10 #min value for plotting 95% bats
# N_min_int <- 10 #min value for plotting boot means
#
# Boot_sums <- Boot_sums %>%
# mutate(low = ifelse(n.obs>=N_min_int, low,NA )) %>%
# mutate(high = ifelse(n.obs>=N_min_int, high, NA )) %>%
# mutate(mean = ifelse(n.obs>=N_min_mean, mean, NA ))
#
#
Faceted bat plots
#code for faceted plots
facet_plot_bats <- function(d=d, ranked= FALSE, Boots = FALSE){
# sample size per species
n <- d %>% group_by(Species_full) %>% summarize(n = length(unique(ID)))
###If Boots= TRUE, calculate bootstrapped values of group
if(Boots){
#bootstrap playback values within each bat spp
Boot_sums <- d %>%
filter(Call_type !="NA") %>%
mutate(effect= paste(Species_full, Playback, sep="_")) %>%
boot_ci2(d = ., y = .$Response, x = .$effect, perms = 5000, bca = F) %>%
mutate(Species_Playback = effect) %>% #keep combined column in new name
separate(effect, into=c('Species_full', 'Playback'), sep="_", remove=F) #re separate combined column
# add n spp
Boot_sums <- left_join(Boot_sums, n) # n per spp.
# delete bootstrapped values below minimum sample size for plotting (could also do by total bat n, not observations per each playback)
N_min_mean <- 10 #min value for plotting boot means
N_min_int <- 10 #min value for plotting bootstrapped 95% intervals
Boot_sums <- Boot_sums %>%
mutate(low = ifelse(n.obs>=N_min_int, low,NA )) %>%
mutate(high = ifelse(n.obs>=N_min_int, high, NA )) %>%
mutate(mean = ifelse(n.obs>=N_min_mean, mean, NA ))
}
#Set factor levels for species from high to low mean response for plotting (doesn't apparently work)
Species_levels <- levels(fct_reorder( d$Species_full, -d$Response, .fun='mean'))
# various failed attempts to make the facet_labels read: "italic(species name ) n= Sample size"
#part of the problem is that I haven't been able to reorder the facets (by changing the factor levels) outside the facet_wrap() command
#appender <- function(string, suffix) paste0(string, suffix)
#appender <- function(string) paste0(string, distinct(subset(Boot_sums, Boot_sums$Species_full == string), n))
# returns all n, not the n of facet
#p + facet_wrap(~am, labeller = as_labeller(appender))
#label facets
lev <- n$Species_full
lab <- paste(n$Species_full, " (n = ", n$n , ")", sep = "" )
names(lab) <- lev
#make dotplot for bats
d %>%
filter(Call_type !="NA") %>%
mutate(Species_full = factor( Species_full, levels = Species_levels)) -> d
#Why don't these reorder facets? infuriating
d$Species_full <- fct_reorder( d$Species_full, -d$Response, .fun='mean')
#reorders the x-axis from high to low mean response if ranked = TRUE
if(ranked){
plot <-
d %>%
mutate(Playback = fct_reorder(Playback, Response, .fun='mean')) %>%
ggplot( aes(x = reorder(Playback, -Response), y = Response))
}else{
plot <- d %>%
ggplot(aes(y=Response, x = Playback))
}
plot <- plot +
facet_wrap(~factor(Species_full,
levels = Species_levels), #order facets by mean response
#facets = ~Species_full,
ncol = 1,
#labeller = as_labeller(appender),
#labeller = as_labeller(lev),
scales = ifelse(ranked, "free_x", "fixed")) +
geom_dotplot( aes(color = Call_type,
fill = Call_type ) ,
binaxis = "y",
stackdir = "center",
method = "histodot" , # fixed bin width dotdensity is default
binwidth =1/40, #max bin width - proportion of range of data
alpha = .7,
stackratio = 0.5, # doesn't change dot size #.5 .9 for seperate good for together plot
dotsize = 6 , #(previously not specified) #2 good for collective plot 8 or 12 good for single spp. plots
# width = .2
# right = TRUE
) +
scale_fill_manual( 'Call type', values = CallCatColors) +
scale_color_manual('Call type', values = CallCatColors) +
# colScale_CallCat + # fill palette
# colScale_color_CallCat + #color palette
coord_cartesian(ylim = c(0, 5)) + #standardize y axis
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_rect(fill = 'white'),
axis.line = element_line(colour = "black"),
text = element_text(color = "black", size = 10),
plot.title = element_text(face = "italic", size = 11),
legend.position = "top",
axis.title.x = element_blank(),
axis.text.x = element_text(angle = 63, hjust=1, face = "italic"),
axis.text = element_text( color = "black", size = 10),
#axis.title = element_text(size = 16),
) +
#dark_theme_noline +
xlab("Playback") +
ylab("Response scores")
#add bootstrapped means if Boots = TRUE
if(Boots){
plot = plot + geom_point(data = Boot_sums,
aes( x= Playback,
y = mean),
color = "black",
shape = 4,
size = 1 ) +
#if more than n individuals, add bootstrapped 95% CI's to plot
geom_errorbar( data = Boot_sums,
aes( x=Playback,
ymin = low,
ymax =high,
y = NULL ),
color = "black",
width=.2,
size=.25)
}
# #make pdfs of plots (85, 129, or 177 mm width to fit in 1-3 columns dimensions)
# ggsave("testplot.jpg", width = 177, height =800, units = "mm")
plot = plot + ggsave(filename = paste("faceted_", ifelse(ranked,"ranked_", "unranked_"), ifelse(Boots,"CIs_", "no_CIs_"), "all_bats_", "responses_to_playbacks", ".jpg", sep = ""), width = 177, height =800, units = "mm" )
#add plot to list
list(plot)
}
unranked_plot <- facet_plot_bats(d = g_u, ranked = FALSE, Boots = FALSE)
## `summarise()` ungrouping output (override with `.groups` argument)
unranked_plot
## [[1]]
#
# ranked_plot <- plot_bats(d= g_u, ranked = TRUE)
# ranked_plot
#TO DO:
#make facet titles italic and add N's (use hum_names strategy below?)
#collapse response scale 3 into 2 (not apparently informative)
#Low priority: ranked isn't working right (uses same order for all facets). Fix?
#pull out dotplot variables into function so can change for each plot as makes sense
#low priority: put Bootsums into function so not hanging out
#put the other plotting variables into the function
#respace dots
####### IGNORE ME. more workspace for figuring out how to add sample sizes to label facets
# #Create an object using as_labeller(). If the column names begin with a number or contain spaces or special characters, don't forget to use back tick marks:
# # Necessary to put RH% into the facet labels
# hum_names <- as_labeller(
# c(`50` = "RH% 50", `60` = "RH% 60",`70` = "RH% 70",
# `80` = "RH% 80",`90` = "RH% 90", `100` = "RH% 100"))
# #Add to the ggplot:
# ggplot(dataframe, aes(x = Temperature.C, y = fit)) +
# geom_line() +
# facet_wrap(~Humidity.RH., nrow = 2, labeller = facet_names)
#make list that reads "Species full = "Species full (n= N)
# facet_names = n %>%
# mutate(facet_names = paste(Species_full, " = " , Species_full, " (n = ", n, ")" , sep = "" )) %>% select(facet_names)
#
#
# # figure out how to add quotes around Species full (n= n_ )
# n<- n %>%
# mutate(facet_names = paste(Species_full, " (n = ", n , ")", sep = "" )) =
#
# # labeller = as_labeller(facet_names)
# labeller = label_bquote(italic(.(sex))~subjects))
#
# labeller = label_bquote(italic(.(Species_full))~subjects))
# # .
# length(unique(.$ID)
# abeller = labeller(cyl =
# c("4" = "condition: 4",
# "6" = "condition: 6",
# "8" = "condition: 8"))
#labeller = label_bquote(italic(.(factor(Species_full, levels = Species_levels)))),
#labeller = as_labeller(facet_names),
# lev <- n$Species_full
# lab <- paste(n$Species_full, " (n = ", n$n , ")", sep = "" )
#
# names(lab) <- lev
#
#
# label_facet <- function(original_var, custom_name){
# lev <- levels(as.factor(original_var))
# lab <- paste0(lev, ": ", custom_name)
# names(lab) <- lev
# return(lab)
# }
#
# facet_wrap(~cyl, labeller = labeller(Species_full = label_facet(original_var = n$Species_full, custom_name = "grouping")))
#
# facet_wrap(~cyl, labeller = labeller(cyl = label_facet(df$cyl, "grouping")))
###### IGNORE ABOVE ################
plot with only top
Individual plots, one bat per plot, ranked and unranked x axis
#Y= response, X= playback, within bat spp. boot_ci()
# store dataframes and plots in lists to be able to make a multiplot or grid object- see how plots look together
# remove all the legends
#try making a multiplot situation (ggarrange from egg package, ?align_plots from the cowplot package?)
#make a version where X isn't reordered every time, so just one X on bottom,
#
#response score- 3 has so few data points- doesn't taper with more- combine 3 and 2?
#clean this up (get rid of excess lists and whatnot)
#pull unique species names
Spp <- unique(g_u$Species)
#initialize lists for plots and boot output
####fixes#####
# pull out ranked from function, test if function works on its own x
#try adding in lists, internally x
#check what output of mutate is
# maybe pull out x and y? response and playback? to make function more universal?
##############
#make seperate plots and bootstrapped 95% cis for each bat
# ranked = X axis ranked from greatest to smallest mean response
# N_min_mean = minimum # for plotting boot means
# N_min_int = minimum # for plotting boot intervals
plot_bats_separate <- function(d=d, ranked= TRUE, N_min_mean=10, N_min_int=10) {
#initialize lists for plots and boot output
myplots <- list()
bootsummary <- list()
#loop through bat spp
for (i in 1:length(Spp)) {
#pull out one bat species at a time
data <- d %>%
filter(Call_type !="NA") %>%
filter(Species ==Spp[i])
# # of bats in group
n <- length(unique(data$ID))
#get bootstrapped means and 95% confidence intervals for the response to each playback in n> 10
if(n >=N_min_mean){
Boot_sums <- boot_ci2(d = data, y = data$Response, x = data$Playback, perms = 5000, bca = F) %>%
rename(Playback = effect)
}
#store bootsums in list
# bootsummary[[i]] <- Boot_sums
#make dotplot for bats
plot <- data %>%
filter(Call_type !="NA")
#reorders the x-axis from high to low mean response if ranked = TRUE
if(ranked){
plot <- plot %>% mutate(Playback = fct_reorder(Playback, Response, .fun='mean')) %>%
ggplot( aes(x = reorder(Playback, -Response), y = Response))
}else{
plot <- plot %>%
ggplot(aes(y=Response, x = Playback))
}
plot <- plot +
geom_dotplot( aes(color = Call_type,
fill = Call_type ) ,
binaxis = "y",
stackdir = "center",
method = "histodot" , # fixed bin width dotdensity is default
binwidth =1/40, #max bin width - proportion of range of data
alpha = .7,
stackratio = 0.9, # doesn't change dot size #.5 good for together plot
dotsize = 8 , #(previously not specified) #2 good for collective plot 12 good for single spp. plots
# width = .2
# right = TRUE
) +
colScale_CallCat + # fill palette
colScale_color_CallCat + #color palette
coord_cartesian(ylim = c(0, 5)) + #standardize y axis
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_rect(fill = 'white'),
axis.line = element_line(colour = "black"),
text = element_text(color = "black", size = 10),
plot.title = element_text(face = "italic", size = 11),
legend.position = "right",
axis.title.x = element_blank(),
axis.text.x = element_text(angle = 63, hjust=1, face = "italic"),
axis.text = element_text( color = "black", size = 10),
#axis.title = element_text(size = 16),
) +
#dark_theme_noline +
xlab("Playback") +
ylab("Response scores") +
ggtitle(paste(unique(data$Species_full)," (n = ", length(unique(data$ID)), ")", sep ="") ) #makes title current group, n is number of bats
if(n >=N_min_mean){
plot = plot + geom_point(data = Boot_sums,
aes( x= Playback,
y = mean),
color = "black",
shape = 4,
size = 1 )
}
#if more than n individuals, add bootstrapped 95% CI's to plot
if(n >=N_min_int){
plot = plot + geom_errorbar( data = Boot_sums,
aes( x=Playback,
ymin = low,
ymax =high,
y = NULL ),
color = "black",
width=.2,
size=.25)
}
# #make pdfs of plots (85, 129, or 177 mm width to fit in 1-3 columns dimensions)
plot = plot + ggsave( filename = paste("individual_", ifelse(ranked,"ranked_", "unranked_"), unique(data$Species_full), "responses_to_playbacks", ".jpg", sep = ""), width = 177, height = 55, units = "mm" )
#add plot to list
myplots[[i]] <- plot
#print(plot)
# print(Boot_sums)
}
names(myplots) <- Spp
myplots
}
#unranked plots
unranked_individ_batplots <- plot_bats_separate(d = g_u, N_min_mean = 10, ranked = FALSE)
#unranked_individ_batplots
#ranked plots
ranked_individ_batplots <- plot_bats_separate(d = g_u, N_min_mean = 10, ranked = TRUE)
ranked_individ_batplots
## $Trin_nice
##
## $Lonc_auri
##
## $Micr_hirs
##
## $Phyl_hast
##
## $Loph_silv
##
## $Glyp_davi
##
## $Gard_cren
##
## $Trac_cirr
##
## $Micr_micr
##
## $Arti_jama
##
## $Tona_saur
##
## $Phyl_disc
##
## $Loph_bras
##
## $Lamp_brach
##
## $Micr_schm
combined plot with ranked x axes
#make plot with aa bats and ranked X axis
ggpubr::ggarrange(plotlist = ranked_individ_batplots[1:15], ncol = 1, common.legend = TRUE)
ggsave("AA_test_plot_allbats_long.jpg", width = 177, units = "mm", height = 800)